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Abstract. We introduce a new micro-macro decomposition of coUisional kinetic 
equations in the specific case of the diffusion Umit, which naturally incorporates 
the incoming boundary conditions. The idea is to write the distribution function 
/ in all its domain as the sum of an equilibrium adapted to the boundary (which 
is not the usual equilibrium associated with /) and a remaining kinetic part. 
This equilibrium is defined such that its incoming velocity moments coincide 
with the incoming velocity moments of the distribution function. A consequence 
of this strategy is that no artificial boundary condition is needed in the micro- 
macro models and the exact boundary condition on / is naturally transposed 
to the macro part of the model. This method provides an 'Asymptotic pre- 
serving' numerical scheme which generates a very good approximation of the 
space boundary values at the diffusive limit, without any mesh refinement in the 
boundary layers. Our numerical results are in very good agreement with the 
exact so-called Chandrasekhar value, which is explicitely known in some simple 
cases. 



1. Introduction 

The development of numerical methods to solve multiscale kinetic equations has 
been the subject of active research in the past years, with applications in various 
fields: plasma physics, rarefied gas dynamics, aerospace engineering, semiconduc- 
tors, radiative transfert, ... The general problem is to construct numerical schemes 
that are able to capture the properties of the various scales in the considered system, 
while the numerical parameters remain as independent as possible of the stiffness 
character of these scales. 

There is a huge literature dealing with the construction of such schemes, and in 
particular it is now a usual challenge to design the so-called Asymptotic Preserving 
(AP) [18j numerical schemes for kinetic equations which are consistent with the 
kinetic model for all positive value of the Knudsen number e, and degenerate into 
consistent schemes with the asymptotic models (compressible Euler and Navier- 
Stokes, diffusion, etc) when e — )• 0. Among a long list of works on this topic, we 
can mention [30l [3ll HSl [20l HSl [23l El for stationary problems and |22l [2ll [21 [271 
[25l|26l[l6l|33l[35l[IIl[5l[l[29]for time dependent problems. 

Here, we are interested in boundary value problems and our aim is to develop 
a strategy to construct numerical schemes that are AP in the usual sense and are 
also able to deal with the space boundary conditions and kinetic boundary layers. 
In this paper, we focus on the diffusive scaling but our strategy could be extended 
to other asymptotics. We emphasize that in general, the known numerical schemes 
induce at the limit e — )• some numerical boundary conditions which are different 
from the theoretical limiting boundary values. In most cases, this does not affect 
only the boundary layer, but also the solution inside the domain. 
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Boundary conditions in the construction of AP schemes for transient kinetic equa- 
tion in the diffusive scahng have been taken into account in several works. For 
instance, in |24| [25| 126) . Klar uses an approximation of the solution to the Milne 
problem which is injected as a boundary condition for the evolution problem. The 
boundary value used by Jin, Pareschi and Toscani in |22| is also derived from the 
asymptotic analysis of the problem as the Knudsen number e goes to zero. Of 
course, taking the boundary condition from the limiting model cannot provide a 
good approximation of the problem for all regime, in particular for large values of 
£. Our goal here is to develop a new approach in which the boundary conditions 
are not extracted from the limiting model but are directly derived from the original 
kinetic boundary condition, at any time. Our numerical results will be compared 
with those obtained with the approaches in |24[ l22l [33] . 

To this purpose, our starting idea is to extend the micro-macro decomposition 
of [21 |33l IS] in order to incorporate in a natural way the exact inflow boundary 
conditions. The general micro-macro decomposition method consists in writing the 
kinetic equation as a system coupling a macroscopic part (say a Maxwellian) and 
a remaining kinetic part. The main advantage of this method is its robustness and 
easy adaptability, since it can be applied to various asymptotics (fluid, diffusion, 
high-field, etc) and to a large class of collision operators. However, the general 
problem of space boundary condition and boundary layers has not been completely 
solved in this approach. When incoming boundary conditions on the distribution 
function are prescribed, it is clear that this cannot be translated into a boundary 
condition on the macro part (say the total density) in an exact way since the dis- 
tribution function is only known for incoming velocities. Such boundary conditions 
for both macro and micro parts of the distribution function are needed for the com- 
putation of fluxes. In this case, an artificial boundary condition of Neumann type is 
used in |33) and leads to some unsatisfactory approximation of the Chandrasekhar 
value at the diffusive limit. 

Our first aim is to develop a new micro-macro decomposition of collisional kinetic 
equations which naturally incorporates the exact space boundary conditions (BC). 
The second task in this work is to show how to use this new formulation at the 
boundary in order to capture the right boundary conditions and boundary layers in 
the limit of a small Knudsen number. 

To achieve this goal we proceed as follows. We consider a kinetic equation with 
a linear collision operator in a diffusive scaling on a bounded domain, subject to 
inflow boundary conditions. Our first idea is to decompose the distribution function 
/ in its domain as the sum of a Maxwellian part adapted to the boundary (which 
is not the usual Maxwellian associated with /) and a remaining kinetic part. This 
Maxwellian is defined on the whole domain such that its 'incoming' velocity moments 
coincide with the 'incoming' velocity moments of the distribution function. In this 
way, the exact boundary condition can be directly injected into the micro-macro 
system. The second idea is to use a suitable approximation of the space derivative 
of the density at the boundary which garantees a balance between the incoming and 
outgoing fluxes at this boundary. Important consequences of this strategy are the 
following. No artificial boundary condition is needed in the micro/macro models 
and the exact boundary conditions on / are naturally shared by the macro part 
and the kinetic part. A further fundamental property of the strategy developed 
here is that it provides AP schemes inside the physical domain and a very good 
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approximation in the space boundary layers. Indeed, at the end of this paper, we 
provide various numerical tests which illustrate this property. We emphasize that 
the numerical boundary value given by our scheme when e is small is very close to 
the theoretical boundary value derived from the Chandrasekhar function, without 
injecting this value in our scheme. The strategy developed in this paper has been 
announced in a preliminary note [32], where we showed that this method can also 
be applied in principle for a hydrodynamic scaling. We will explore in the future 
several extensions of this work, in particular the possibility of dealing with more 
complex collision operators. Our approach shares some similarities with the method 
developed in |10| [12] where a partial moment system with entropy closure was used 
to approximate the original kinetic equation. 

The paper is organized as follows. In Section [2] the general strategy is described 
at the continuous level. We introduce a new micro-macro decomposition which 
is adapted to boundary value problems for kinetic equations. In Section [3j we 
use this decomposition to design a numerical scheme which is AP in the diffusion 
limit. In a first step (Subsection 3.1), we show how this formulation allows to 
inject directly the inflow boundary condition in the numerical method, avoiding 
in this way any artificial boundary condition. In a second step (Subsection 3.2), 
we introduce a specific (still consistent) discretization of the spatial derivative of 
the density in order to correctly capture the boundary condition and the boundary 



layer. The properties of our scheme are summarized in Proposition 3.2 In Section 



|4j we provide various numerical tests to validate our approach. We test different 
boundary conditions (with or without boundary layers) and different linear collision 
operator (with local and non local collision kernels). Finally, we end the paper by 
a conclusion and some perspectives. 



2. A boundary matching micro-macro decomposition 

2.1. The diffusion hmit of kinetic equations in bounded domains. Let Q. 

be a domain in the position space and let V he a. domain in the velocity space 
M"^, V being endowed with a measure d^. At any point x G dVl of the boundary of 
il, we denote by n{x) the unitary outgoing normal vector to dO.. We then consider 
the following linear transport equation written in a diffusive scaling: 

edtf + =^^Lf, t>0, {x,v)£VLxV, (2.1) 

with initial condition 

f{^,x,v) = finit{x,v), {x,v)(^VLxV (2.2) 

and with inflow boundary conditions 

f{t,x,v) = fb{t,x,v), t>0, {x,v) e dfl X V such that V ■ n{x) < 0. (2.3) 

The unknown / is the distribution function of the particles that depends on time 
t > 0, on position x S and on velocity v € V. 

In (2.1 ), the linear collision operator L acts on the velocity dependence of / and 
describes the interactions of particles with the medium. This operator relaxes the 
system of particles to an equilibrium £{v), which is supposed to be a positive and 
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even function. For all distribution function h, we shall denote 

, , h(v)dv 
Jy£{v)dv 

Note in particular that {£)y = 1 and {v£)y = 0. We assume that the collision 
operator L is non-positive and self-adjoint in L^{V,£^^dfi), with null space and 
range given by 

Mil) = Span{£:} 

and 

7^(L) = (AA(L))^ = {/ such that {f)y = 0}. 

Hence, an important property of L is the fact that it is invertible on TZ{L), we shall 
denote its pseudo-inverse by L~^. Additionally, we assume the invariance of L under 
orthogonal transformations of M'^. The parameter e > measures a dimensionless 
mean free path of the particles, or also the inverse of the dimensionless observation 
time. 



When e goes to in (2.1), / goes formally to an equilibrium state /o(t,x,f) = 
po{t,x)£{v). The diffusion limit is the equation satisfied by the limiting density po- 
This equation is a diffusion equation and has been obtained in various situations 
thanks to asymptotic expansion in e (Hilbert or Chapman-Enskog expansion), see 
e.g. [SSlElIIliaEi. We have 

dtPo-V^-{KV^po) = with K = -{vL-\v£))^>0, (2.5) 

with a Dirichlet boundary condition 

po{t,x) = pbit,x), t>0, x£dn. (2.6) 

The limit value pi, is usually obtained by a boundary layer analysis and its com- 
putation involves the resolution of a kinetic half-space problem, the Milne problem 
associated with the collision operator L. For x G dQ, let X^'^iv^'^) be the bounded 
solution of the half-space problem 

- V ■ n{x)dyx^''' = Lx^'"", y>0, v £ V, 

X*'"(0,t;) = fb{t,x,v), vn{x)<0, (2.7) 

in which t and x are parameters. Then, under reasonnable conditions of /;,, this 
problem is well-posed |3l [U [TH |36] and we have 

hm x''^y,v)=pb{t,x)£{v). (2.8) 

J/— >+oo 

Notice that, in the special case where the incoming data is an equilibrium state, i.e. 
when fi,{t,x,v) = pb{t, x)£{v), one has X^'^iu^'v) = Pbit, x)£{v) and is independent 
of y. In this case, there is no boundary layer near di} as e — t- 0. For general incoming 
distribution functions, the question of determining pt, requires the resolution of the 



half-space problem (2.7), which may be computationaly demanding. Our aim here is 



to construct a numerical method which gives a good approximation of the boundary 



layer without solving the half-space problem (2.7). 



In order to check the efficiency of our approach, a validation test will be the 
following. We will compare our numerical results with the exact Dirichlet value pi, 



defined by (2.8 ) in the diffusion limit e — )• 0, in special cases where this value is well- 



known. The most simple case is the one-group transport equation in slab geometry. 
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which is such that d = 1, Q = [0,1], the velocity set is V = [—1, 1], dfi = ^dv and 
Lf = {f)v - /■ Then pb is given by 

Pb{t,0)= [ Ko{v)fb{t,0,v)dv, Pb{t,l)= [ Ko{-v)fb{t,l,v)dv, (2.9) 



where 

Ko{v) = ^vH{v), (2.10) 

the function H being the so-called Chandrasekhar function [TJ [H E], see also [6], 
which can be computed numerically from the following formula: 

' 2 Jo v + w 

A good approximation of the function Kq{x) is usually given by 

Ko{v)c^Ki{v) = ^v^ + v, (2.12) 

see 1301 for instance. 



2.2. A new micro-macro decomposition. Let us first recall the micro-macro 
decomposition introduced in |33j and |34| . We decompose 

f = p£ + gi, (2.13) 

with 

P{t, x) = {f)y and gi = f - p£ 
{gi is not necessarily small). We denote by 11 the orthogonal projector in L'^{£^^dp) 
onto the nullspace of L, i.e. 

U<P={<P)y£. 



Inserting the decomposition (2.13) into the kinetic equation (2.1) and applying 11 



and / — n successively (/ being the identity operator), one gets 
dtp+ ^"^x ■ {vgi)v = 0, 

dtgi + ^(/ - n)(^; • V^gi) = ^ [Lgi - eSv ■ V^/))] . 

We now emphasize that in this formulation, the space boundary condition on p 
and gi are not known because they cannot be inferred from the boundary condi- 



tions on / in general. Indeed, (2.3) only gives the function / at the boundary for 
incoming velocities, i.e. such that v ■ n{x) < 0, and it is therefore clear that the 
values of p{t,x) = {f)y cannot be determined a priori on the boundary dQ with 



the only knowledge of fb- Consequently, the micro-macro decomposition (2.13) ap- 
pears as a method for capturing the diffusion limit inside the domain only. The 
numerical scheme constructed in [33] is Asymptotic Preserving in a strict sense 
(i.e. including the boundary) only for incoming distribution data at the equilib- 
rium fb{t,x,v) = pb{t,x)£{v), when there is no boundary layer near dVt. Indeed, 
only a rough approximation of the boundary layer is obtained for non equilibrium 
boundary conditions. 

Our aim in this work is to develop a new micro-macro decomposition which is 
able to incorporate the boundary conditions in an exact way. This will allow us in 
the next section to construct a numerical method with natural boundary conditions. 
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that approximate very well the diffusion limit even for incoming data that are not 
an equilibrium state. 

To this purpose, let us introduce a few further notations. We consider a function 
v) which extends the function n{x) ■ v (defined on dQ x V) to the whole domain 
0,xV. Let us give explicit examples of u}{x,v) for specific geometries. For the unit 
ball centered at the origin, one can take u}{x,v) = x ■ v. For a half plane xi > 0, 
one can choose u}{x, v) = {—vi,0, 0). In dimension one, for the interval [0, 1], one 
can take uj{x, v) = {2x — l)v. Then, for all x £ we split the velocity space V into 
two parts according to the sign of this function uj{x, v): 

V-{x) = {vGV, uj{x,v) < 0}, V+{x) = V\V^{x). 

For all function h{v), we will denote 



(2.14) 



Jy_ h{v)dfi 
Jy_£{v)dfi' 



Iv+ Hv)df^ 
jy^£{v)dix 



and 



£. 



Ilv.h={h)y_ 

The idea is now the following. Instead of looking for an equation on p as usual, we 
seek an equation on the following 'boundary matching' density: 

= (/(i,x,-)>v„ (2.15) 
and perform the corresponding micro-macro decomposition: 

f = pE + g = I[v_f + g. (2.16) 

In particular, we have {g)y_ = 0. When e — )• 0, we know that the solution / of (2.1 ) 
is, at least formally, close to p£ except in initial or boundary layers. Therefore, since 
we have {p£)y = p, the moment p will also be close to p for small e. This shows 
that the new decomposition (2.16) is still a decomposition of / into an asymptotic 
part (macro part) and a kinetic part (micro part). In order to derive the system of 

(2.17) 



equations satisfied by p and g, we first integrate (2.1) on and get 

72 (^9)v. ■ 



1 



dtp + - {v£)y_ 



1 



• ^xP+ - {v ■ Vxg)v_ 



e 



Substracting ( |2.17 ) from equation (2.1) on /, we obtain the equation on g: 

dtg + -{I-IivMv y.g) + - (/ - ny_ ) {v£) ■ V^p =\{I-Uv_) (Lg) . (2.18) 

Finally, we remark that the system ( |2.17 ), (2.18) can be replaced by a more con- 
venient (and still equivalent) system in terms of p = {f)y and g = f — p£ as 
follows: 



dtp + ^ ■ ^xg)v 
dt9+l{I 



0, 



1 



ny_ ) {v ■ V^g) + -{I-Uv_) {v£) ■ V^p 



(2.19) 

;{I-Uv_)iLg). (2.20) 



Note that p is linked to p and g by the relations 

P = P-{9)v^ f = p£ + 9 



{g)y£. (2.21) 



One main interest of this new micro-macro formulation (2.16) is the following: the 
fluxes involved in (2.19) and (2.20) only concern the quantities g or ~p, and not p. 
This means that numerical schemes of this formulation would only need the values 
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of p and g at the space boundary which are completely known from the original 



boundary condition (2.3) on /: 



pit,x,v) = {fb{t,x,-))v_ > 

g{t, X, v) = fb{t, X, v) - {fb{t, X, •))y_ £{v), 



t>0, {x,v)£dn, (2.22) 
t > 0, {x,v) edUxV-. (2.23) 



Note that a similar decomposition as (2.16) can be performed for a fluid scaling as 
well, see |32j . 



3. The numerical method 



In this section, we construct in dimension d = 1 a numerical scheme for (2.19) 



(2.20), see Proposition 3.2 which is uniformly stable in the limit e — )• and provides 



a good approximation of the boundary layers, without introducing any artificial 
boundary condition. To this purpose, we proceed in two steps. First, we present 



in subsection 3.1 a preliminary numerical scheme, see Proposition |3.1[ which has 
the desired uniform stability and does not use any artificial boundary condition. 
However, it turns out that this scheme does not give a good approximation of the 
limiting boundary condition. To remedy this problem, we then introduce our nu- 



merical scheme in subsection 3.2 where a specific discretization is done near the 
boundary. We emphasize that this scheme is able to reproduce an accurate approx- 
imation of the boundary layer without any mesh refinement. 



Let us introduce a few notations. The domain in space is the interval [0, 1]. The 
incoming boundary conditions are 



f{t,0,v) = Mv), for veV.{0) = {v€V,v> 0}, 
f{t,l,v) = fr{v), for v£V.{l) = {v(^V,v< 0}, 



(3.1) 
(3.2) 



the data fi and fr being independent of time here for simplicity. The function 
u{x, v) is a smooth function satisfying 



uj{0,v) 



u;{l,v) 



V. 



(3.3) 



-jyi-j- be a uniform space step. We set 
iAx for i = 0, . . . , N + 1 and Xi_^_i/2 = 



Let At be a constant time step and Ax 
tn = nAt and consider staggered grids X; 

{i + 1/2) Ax for i = 0,. . . ,N. We will construct approximations of the densities 
pf ~ p{tn,Xi) and ~ p{tn,Xi) on the grid {xj}o<j<Ar+i and approximations of 
the remainder 5fJY^/2(^) - 9{in,Xi^ii2,v) on the grid {xj+i/2}o<i<Ar- The velocity 
variable v is kept continuous in this section, and its discretization will be made 
precise later on. At the time t = 0, the initial condition (2.2) induces naturally the 
initialization 



Pi — {finit{Xi, ■))y , Pi — {finit{Xi, •)) 



V- ' 



9i+l/2 



{I -UY-)finit{Xi 



-1/2; ' 



(3.4) 

Note that, by construction, the function p is known at the boundary x = 0, x = 1 



and can be deduced from the boundary conditions (3.1), (3.2): 



PQ 



/y_(o)'^(^)^/^' 



Pn+1 



Iv-il) fr{v)dp, 



(3.5) 
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At the boundary, the function g is also known for entering velocities. We introduce 
two additional unknowns g^iv) and g'^_^_l{v) for v €z V. The incoming boundary 
conditions yield 



9^{v) = 



= fr{v)-p^,,£ iorveV-il). 



(3.6) 



For outgoing velocities, these functions will be provided by our numerical scheme, 
when their values are required, and no artificial condition will be needed. This is 



explained in Subsection 3.2 



3.1. A first attempt without specific treatment at tiie boundary. We as- 
sume that we know /o" and pf for i = 1, . . . ,N, and ^ ~ 0, . . . A^. We 
proceed by recursion and explain how to compute these same quantities at time 
tn+i- First, for i = 0, - ■ ■ ,N, we discretize (2.20) by a semi-implicit upwind scheme: 

1 , , _ .1,9., 



n+l 
9i+l/2 



9 



i+1/2 



At 



+ -(/-nyj 



V 



i+1/2 



+ -{I-UvJ{v£) 



Ax 

Pi+1 Pi 

Ax 



9i-l/2 _5i+3/2 

- — \-v — 



9. 



i+1/2 



1 



Ax 



ny_)(L5"+i 



i+1/2' 



(3.7) 



where = max(f,0) and v = min(?;,0). Note that the source term involving 
dxP is discretized by a centered finite difference scheme. We observe that g^i^2 ^^'^ 
9n+3/2 incoming velocities can be extracted from the relations 



9o 



9-1/2 + 9i/2 



9 N+l 



9n+1/2 + 9n+3/2 



(3.8) 



2 ' 2 
since g^^ ^i^'^ 5]v+i/2 known and since and (7]y+i known for incoming 
velocities by (3.6). Note also that ~Pq and pjv+i known by (3.5). Therefore, the 
scheme (3.7) enables to compute all the desired values g^^y2 ^ ~ 0, . . . A^ (see 
below in the proof of Proposition 3.1). 

It remains to compute p^ and p^ for i = 1, . . . , A^. To that purpose, we consider 
the following discretization of (2.19). For z = 1, . . . , A^, we set 



^n+l 



At 



+ 



1 



n+l 
9i+l/2 



n+l ' 
9i-l/2 



Ax 



0. 



V 



These relations clearly provide all the desired values for i 
compute all the values according to (2.21), by setting 



Tzn+l 



p: 



n+l 



9i+l/2 ^ 9i 



-1/2 



(3.9) 

, A^. Then we 
(3.10) 



V 



for i = 1, . . . , A^. The recursion is now complete and can be iterated. Observe that, 
in this first attempt, we do not use the outgoing values of g^ and fi'Jv+i ™ order to 
compute p^~^^ and ^^^^2 ™side the domain. 

We summarize in the following proposition some properties of this scheme. 



Proposition 3.1 (Formal). The numerical scheme (3.7), (3.8), (3.9), (3.10) with 
initial and boundary conditions (3.4), (3.5), (3.6) determines all the values p^ and 
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for i = 1, . . . , iV, and 5-[Yi/2 /^^^ ^ = 0, . . . A^, 
scheme is stable under the CFL condition 



at any time interation n. This 



At<C (Ax^ + eAx) , 
where C is a constant independent of e, At and Ax. For all fixed e > 0, this scheme 



is consistent with the initial and boundary value problem (2.1), (2.2), (2.3). More- 
over, for all fixed At and Ax, the scheme degenerates as e — )• into a discretization 
of the diffusion equation (2.5). More precisely, for i = 1, 
0, and is defined by the scheme 

2p« + pt 



, A^, we have pf ^ 



as e 



p? 



~n 



At 



Ax2 



0, for i = l,...,N, 



where k > is given in (2.5) and where the boundary values are 
{vL-\l 



Po 



PN+1 



{vL~^I-U){v-fr))v 
{vL-^{v£))v 



we 



miv-'fe))v 
{vL-\v£))v 
for all n. We recall that Ilk = {h)v£- 

In the particular case where V = [—1, 1], £ = 1, dp = ^dv and L = II — I 
obtain k = ^ and 

Po= K2{v)fe{v)dv, pN+i -- 
Jo 

with 

K2{v) = 3v'^. 

For instance, if /^(v) = v, we get the value po 
limit ing d iffusion equation, whereas the right value can be computed from (2.9) 

and (2.10): pb{0) = ^ j^v^H{v)dv ~ 0.7104. The value obtained by this scheme 
should clearly be improved. This is the goal of the next section. Note that the 
same boundary value 0.75 is obtained by the schemes proposed in |33| and |22| . 
which both involve some artificial boundary conditions. We emphasize that, by 



K2{-v)fr{v)dv. 

(3.11) 

0.75 as boundary value for the 



construction, our scheme only involves the natural boundary conditions (3.5) and 



(3.6) derived from fi and fj. and does not need any artificial condition. 



Proof of Proposition 3.1 

We first check that all the values are uniquely defined by the scheme. The algorithm 



has been described above and it only remains to show that (3.7) allows to compute 
9^+1/2 ™ terms of quantities at time n. To this aim, we simply remark that the 
calculation of §^^-^^2 done by solving a linear system of the form 



At 



Since the operator 



I-L]g = h, with L = (I - UvjL{I - 



(3.12) 



-L and 



2 

then the operator ^/ 



L is self-adjoint and nonnegative, so is the operator 
L is invertible, for e > 0. 
For the stability CFL condition, a similar proof as in \33\ |35] can be done. We 
only sketch the main arguments of this proof and refer to these works for the details. 
Roughly, when e = 0(1), the convection term in the equation (3.7) prevails and 



induces the standard stability condition At < CeAx. When e <^ 1, the system 
behaves as the diffusion limiting equation (see below) and this results in the usual 
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stability condition At < CAx^. In the general case, the CFL condition is expressed 
as a combination of these two hmiting situations. 

Now, we analyze formally the scheme when e — )• 0, the parameters At, Ax 
being fixed. First, we remark that the operator (/ — Ilv_)L is invertible from 
{/ • {f)v- = 0} to itself. More precisely, one has 

(/ - Uv_)Lg = h, with ge{f: (/)y_ = 0}, 

if, and only if, 

g = {I- Ilv_)L-\l - U)h, with hG{f: = 0}. 

0{e), for i 



Now, by induction, it is clear from (3.7) that g^^i^2 
and n > 1, and then, from (3.9), we deduce that /o" = 0(1) and 
i = 1, . . . , N and n > 1. Reinserting these estimates in (3.7), we get 

Pi+i 



n+l 



e{I-Uv_)L-\v£) 

for i = 1, . . . , — 1, and 



p7 



Ax 



+ O(e') 



0,...,A^ 

O(l) for 



(3.13) 



Pi 



Po 



Ax 



Ax 



(I-U) 



= 5(i-nv._)L-i 



.Sr^^)^liI-U)i.-g^^,,) 



Inserting these equations into (3.9) gives, after some computations. 



pr' - p'l 



At 



Ax2 



0{e), for i = 2,...,iV-l 



and 



At Ax2 ^ 



Ptv"^ ~ p1<! Pn+ 1 ^Pn + Pn-1 

Ki 



At 



Ax2 



Oie), 



the quantities and p'^^i being defined in Proposition 3.1 



□ 



3.2. The numerical scheme. In this section, we construct a numerical scheme 
which is Asymptotic Preserving inside the domain and which provides a good ap- 
proximation of the exact value at the boundary. To this purpose, we proceed as 
follows. Inside the space domain, we use the same numerical scheme as the one de- 
scribed in the previous section. At the boundary, we propose a specific treatment, 
without using neither artificial boundary condition nor any mesh refinement. 

Our approach is based on the following idea. We consider the spatial derivative 
of the density at the boundary as an additional unknown in the scheme. It is known 
indeed that the density becomes stiff (its derivative is of order ^ in the boundary 
layer) and the analysis of the boundary layer requires usually a rescaling process. 
Here, our strategy enables to avoid such rescaling. In the boundary cells, we simply 
determine the slope of p by ensuring the balance between the incoming and outgoing 
half- fluxes of /. We recall indeed that, if x is a solution of the half-space problem 
(2.7), then one has dy{vx)v = and then, from (2.8), we deduce that {vx)v = 

vxdp = 



and then, at y = 0, we have 
the time evolution context. 



V+(0) 



Iv (0) ^fb^P- shall use this idea in 
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Let us now describe our numerical scheme. We assume that we know and 
for i = 0, . . . , N + 1, and g^j^i^2 i = 0, ■ ■ ■ N . We also assume that we know the 
approximation of the function g at the boundary, i.e. Qq and g^_^_l for all v £ V. 
We will now show how to compute all these quantities at time tn+i- 



First step: determination of the quantities at time iteration n + 1, near the bound- 
ar^es: <+S and g-/_^\^„ g^+i 



First, we use a centered scheme for the calculation of g^J^ and g^lf^i/2' 

, l„_n..,f„?I^) ,3,14) 

n+l\ 



At e' '-'V Ax 



5^+1/2 ^ ^(j-nv^jf /^+i-g^ ') (3.15) 



At 'V Ax 



-'■/'r TT \f„.c\ f Pn+1 - Pn\ _ ^ (T TT 



+ -{I-Yiy_){v£) ( ^^^^^^^ = -Ilv.){Lgl^y,> 



where g^ and (7]^ are determined by interpolation 



9i/2 + 53/2 „ _ 9n -1/2 + 9n+i/2 

2 ' 2 • 

It remains to determine gQ~^^ and 5]\r!^\- incoming velocities, these functions 



are prescribed thanks to (3.6): 



lF_(0)5r'W = ly_(o)(/K^)-Po'?), (3.16) 
ly_(i)5]^X\W = Iv.(i) if r{v)-p^+,£), (3.17) 
where ly^ is the characteristic function of a set A. For outgoing velocities, we dis- 



cretize the transport equation (2.20 ) in a suitable way. The specific point here is that 
we consider the derivative of p at the boundary as an unknown. We approximate 
£^^(0) by the quantity A" and £^^{L) by the quantity A", which are computed 
below. For v < only, we write 

+ -il - ny.) (v ''^'~f'^' ] = -M - ny_)(L<+i) (3.18) 



and, for f > only, we write 
At 



+ }.^J _ n._) r/jli^llZlE^ \ J_„ _ n,,_)(L«\)(3.i£.) 
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In these equations, we have used the quantities g^Yji ^^"^ 5']v+3/2 '^^ich are again 
defined by extrapolation: 



5-1/2 — ^^O ~ 51/2 ^^I'i 57V+3/2 — "^Qn+I ~ 9% +1/2 ■ 



(3.20) 



Let us give a sense to these two equations (3.18) and (3.19) which are both written 
on V+. We only deal with the case of (3.18), the other equation can be treated 
similarly. This equation takes the form 

^/-L^5^+l = ly+/i, With L= (/-ny_)L(I-ny_). 

We recall that Qq'^^ is already known on V- by (3.16). Hence, we have to solve the 
following linear system on the unknown '^v+Qq 

-2 



n+1 



where the right-hand side is given. To see that this system is well posed, it suffices 
to recall that the operator L being self-adjoint and nonpositive, the bilinear form 
associated to — L is symmetric positive definite on the set of functions defined on 
V. Hence this bilinear form is also symmetric positive definite on the set of functions 
defined on V+, and consequently the operator ly^ — l\Tlv^ is invertible on this 
set: Eq. (3.18) admits a unique solution ly^^Q^"^. From (3.18), we get 

^1 



n+1 



At 



90 



X^+Hl-Uv_){v£) 



-e{I-Uv_){v 



9l/2 5_i/2 

Ax 



n+l 



(3.21) 



In order to determine A""*""^ and A""*""^, we discretize (2.19) at x = and at x = 1: 

Po 



Po^' 



At 



, 1 (J^ 



^-1/2 \ 



Ax 



/ 



0, 



V 



Pp+i Pn+1 I 1 
At e\ 



n+1 _ n+1 ' 
yN+3/2 yN+1/2 

Ax 



0. 



Using (2.21), (3.5) and (3.20), we deduce the following two equations: 

9o" 



at' 



At 



9 N+1 9 N+1 



V 



, 2 / ^r/V 

^eV^A^ 



9^' 







V 



At 



2 / ^AT+l 5]^+\/2\ 



0. 



(3.22) 



(3.23) 



Inserting (3.21), together with (3.16 



A""*"^, from which we get g^'^^ using (3.21). Finally, we compute Pq"*"^ using (2.21) 



^n+l 



into (3.22) enables to compute explicitely 



ji+i 



and (3.5). We proceed similarly with (3.19) and (3.23) to compute A""*"^, 5]v+i ^^"^ 



n+1 

Pn+1- 
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Second step: determination of the quantities at time iteration n + 1, inside the 
domain. 



In this second step, we focus on the interior of the domain and compute the values 
p^+\ p^+i for i = 1, . . . , iV and g^^^^^ for i = 1, . . . , - 1. 
For i = 1, . . . , — 1, we write 



n+l 



i+l/2 5'j+l/2 
~^t 



Then, for i = 1, 



+ 



n 

+ y i+l/2 



5" 



+ 




'i-l/2 , _5i+3/2 
+ V — 



9. 



i+l/2 



Ax 



i+l/2l 



, N ^ we write 

.n+l _n 



At 



+ 



n+l 

yi+1/2 



n+l ' 

yj-1/2 



Ax 



0. 



(3.24) 



(3.25) 



Note that the values g^tc^ and /9 have been computed in the first step. 



Proposition 3.2 (Formal). The numerical scheme (3.14), (3.15), (3.18), (3.19), 
(3.22), (3.23), ( 3.24[ ), ( |3.25 ) with initial and boundary conditions (3.4), (3.5), (3.6) 
determines all the values and for i = 0, . . . , + 1, and i?|!)_i/2 ■Z^"'" ^ = 0, . . . A^, 
at any time interation n. This scheme is stable under the CFL condition 

At<C (Ax2 + eAx) , 

where C is a constant independent of e, At and Ax. For all fixed e > 0, this scheme 



is consistent with the initial and boundary value problem {2.1), (2.2), (2.3). More- 
over, for all fixed At and Ax, the scheme degenerates as e — t- into a discretization 
of the diffusion equation (2.5). More precisely, for i = 1, 
as e 



, A^, we have ^ 



0, and is defined by the scheme 
Pi — Pi Pi+i "^Pi ^ Pi-i 



At 



Ax2 



0, for i = l,...,N, 



where k > is given in (|2.5|) and where the limiting boundary values take the form 

KL{-v)fr{v)dv, 



Pq 



V-{0) 



KL{v)ffXv)dv, 



Pn+1 



v-(i) 



(3.26) 



where Kl is an explicit function which depends only on the operator L. The expres- 
sion of this kernel can be deduced from ( 3.32[ ) and (3.33) below. 

In the particular case where V = [—1, 1], £ = 1, d/i = ^dv and L = 11 — /, we 
obtain k = I and we can compute the following expressions for the kernels: 



Po 



Ks{v)fe{v)dv, p^+i 



K'i{-v)fr{v)dv. 



with 



3 2 15 

—V H V 

2 14 



1 

28' 



(3.27) 



For instance, if fi{v) = f, we get the value po = 0.7143 as boundary value for 
the limiting diffusion equation, where we recall that the right value is given by 
Pb{^) = Jq v'^H{v)dv ~ 0.7104. Therefore, this scheme provides a much better 
approximation of the limiting boundary value than the one obtained from the first 
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scheme presented in Sec tion |3.1| On Figure [T] we plotted four curves: the exact 
kernel Kq^v) defined by (2.10), its usual approximation Ki{v) defined by (2.12), its 



approximation K2{v) defined by (3.11) and given by the first scheme presented in 



Section 3.1 and its approximation i^3(u) defined by (3.27). We observe that our 



kernel Ks{v), as well as Ki{v), fit very well with the exact Chandrasekhar kernel 
Kq{v), whereas the kernel K2{v) does not provide a good approximation. 




Figure 1. Exact kernel Kq{v) given by (|2.10) and its approxima- 



tions Ki, K2, K3 defined respectively by (2.12 



(3.11), (3.27) 



Proof of Proposition 3.2 The proof is similar to that of Proposition |3.1[ except that 
we have to deal with the boundary terms. We focus on the left boundary x = 0, the 
right boundary can be treated similarly. Let e > be fixed. Our numerical scheme 
is consistent with the continuous model as soon as the following property holds: 



e 



Pi - Po 
Ax 



+ 0{At,Ax). 



(3.28) 



To show this property, we consider a modified scheme consisting in the same equa- 



tions (|3.14|), ([3.151), (|3.24|), (|3.25|), and where (|3.18|) is replaced by 
9E 



at' 



At 



+ 



1 /n'^+l - 

-(/-nv-_)("f) ^ 



1 



+ -(/-Hv 



9l/2 5^-1/2 

' Ax 



1 



(/-ny_)(L5, 



(3.29) 



for u < 0. This equation allows to compute gQ^'^iv) for u < 0. Since it is clearly 
a consistent discretization of (2.20), the whole modified scheme will be consistent 
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with the continuous model and therefore satisfies 



% 9o_ 

At 



2 / 51/2' - 9o 

+ " \ ^ X — 

e \ Ax 



0{At,Ax). 



V 



Now, let us recall that, in our numerical scheme, the quantity A" was computed 
thanks to (3.18) and (3.22). This enables to conclude that A""*"^ is a discretization 

of and (3.28) holds. 

Let us now analyze the limit e — )• of our scheme near the left boundary. From 



(3.14), we get 



g^l-;=e{I-Uv_)L-' 



v£ 



Pi - Po 
Ax 



Ax 



U-n)Ho") 



+ 0{e'^). (3.30) 



Now, we need to pass to the limit in (3.21). For this, we recall that the bilinear 
form associated to —L = —(I — ny_)L(/ — Xly^) is symmetric and positive definite 
on the set {/ : = 0}. Hence, since the set of functions which vanish on V- 

is included in this set, the operator ly^L ly^ is invertible from the set of functions 
defined on V+ to itself. Thus, from (3.16) and (3.21), we obtain 

iv-9o{v) = Iv- (Mv) - Po£) , 

-1 



X^{I-UvJiv£)-Llv.g^ +0{e) 



(3.31) 



It remains to determine A". From (3.22), we get the equality of the half-fluxes at 
the boundary: 



{vg^)y = 0{e). 



Inserting (3.31) into this equation yields A" = Af + 0{e)^ where 

(-v+{h - Po£) + v-{lv+LlvX' 
Xf = ^ 



Lilv_{h-PoS)) 



-(V+LVj-i [{I-UvJ{v£)] 



(3.32) 



V 



Next, equation (3.25) on p" becomes 
Pi^'-Pl 



At 



Ax^ 



Ax2 



0(e) 



to 



where we used the expression (3.13) for g'^J^ (3.30) for g^^ ■ It remains 

replace g^ in this formula by its expression (3.31), in which A" is given by (3.32). 
We finally get 



Pi 



n+l 



where p^ is given by 
Po = Po-Z 



At Ax2 

vL-Hl-U) [v+{fe-po£)])y 



0{e), 



(3.33) 



1 



vL-Hl-U) 



Xe{I - Uv_ ) (vS) - Llv_ ifi - Po£) 



Finally, we remark that inserting the expression (3.32) of A^ into (3.33) yields the 
existence of the kernel Kl such that (3.26) holds. The proof of Proposition 3.2 is 
complete. □ 
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Let us compute the approximate boundary value given by our scheme in a special 



case of interest, which is tested numerically below. We set V 
and take the collision operator under the form 

'1 



-1,1], d/i 



Lfiv) 



a{v,w){f{w) — f{v))dw with a{v,w) =p{v)p{w), 



1 



p being a given even function. Then straightforward but lengthy computations from 



(3^ and yield 



Po 



K^P{v)h{v)dv, 







where the kernel is given by 



K'^PP{v) = a— + hv + cp{v) + d. 
p(v) 



(3.34) 



The coefficients a, 5, c and d are given by 
1 , 7 + 2a(5 



a = 1. — , 
where we have set 



27^ 



K,a{l + 4q!7) 



Ka(l + 4a7) ' 



d=l 



a 



p{v)dv, (3 



Piv. 



dv, 7 



Piv, 



dv, 



{p{v)y 



dv. 



Notice that in the simplified case where L = 11 — /, i.e. p{v) = 1, then one has 



15 
14' 



h and d = 0, which covers the expression (3.27). 



28 



4. Numerical tests 

We present here some numerical experiments to validate the approach. In our 
tests, the space domain is the interval = [0, 1], the velocity domain isV= [—1, 1] 
and djj, = ^dv. All the integrals in velocity which are involved in the model are 
computed using a Gauss quadrature method. The collision operator will be of the 
following form: 

1 /-i 

Lf{v) = - J ^a{v,w){f{w) - f{v))dw, 

where cr is a given symmetric and nonnegative kernel. We have £{v) = 1. The initial 
data will be zero: finit{x,v) = 0. The function uj{x,v) used in our simulations is 
oj{x,v) = {2x — l)v — x{l — x). In the tests using our numerical scheme, the time 
step is linked to the space grid step as follows: At = — h eAx). 

We compare our numerical results, which will be referred to as LMe, with those 
given by the following schemes: 



The time explicit upwind scheme for the original kinetic equation (2.1), where 
At and Ax are linked to e by standard CFL stability condition. When e > 
10~^, this explicit scheme will be highly resolved to serve as a reference for our 
comparisons. 



ii. The direct discretization of the diffusion equation (2.5) by a standard explicit 



scheme, in the cases where the right boundary condition is known (see e.g. the 



Chandrasekhar value (2.9)). Our results will be compared with those given by 



this scheme when e is very small. 
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iii. The micro-macro scheme of Lemou and Mieussens given in [33j . which will be 
referred to as LMi. 

iv. The scheme by Klar [24j [25], which will be referred to as K. 

V. The scheme by Jin, Pareschi and Toscani ^22], which will be referred to as JPT. 
In all the figures, we plot the density p as a function of space, at different times. 

4.1. Boundary conditions at thermal equilibrium. In this subsection, we 
present some numerical tests in a situation where the incoming boundary condi- 
tion is proportional to the equilibrium £{v) = 1. We recall that, in this case, there 
is no boundary layer in the diffusive limit. 

Example 1. Kinetic regime with isotropic boundary conditions: 

a{v,w) = l, h{v) = l, fr{v) = 0, £=l. 

On Figure [4j we plot the results obtained with our scheme LMe and the reference 
obtained by the explicit scheme at times t = 0.2, t = 0.5, t = 1, t = 2 and t = 4. 
Our scheme LMe is simulated with 50 and 200 grid points in space, whereas the 
explicit scheme is simulated with 2000 grid points. We clearly see that our results 
are in good agreement with the reference simulation in this kinetic regime. 

Example 2. Diffusive regime with isotropic boundary conditions: 

aiv,w) = l, feiv) = l, friv) = 0, £ = 10"^ 

On Figure [5| we plot the results obtained with our scheme LMe and the reference 
obtained by the diffusion equation at times t = 0.01, t = 0.1 and t = 0.4. Our 
scheme LMe is simulated with 50 and 200 grid points in space, whereas the diffusion 
limit is simulated with 200 grid points. We also see here that our results are in total 
agreement with the diffusion limit. Of course, the AP property of our scheme allows 
us to keep At and Ax independent of e. 

4.2. Non equilibrium boundary conditions and boundary layers. In this 
subsection, we present some numerical tests where the boundary data are not at 
equilibrium, which induce a boundary layer in the diffusive regimes. Here, we still 
consider the situation where the cross section a{v, w) is constant, we differ to the 
next section the more general case. 

Example 3. Kinetic regime with non isotropic boundary conditions: 

a{v,w) = l, fi{v) = v, fr{v) = 0, e = l. 

On Figure [6j we plot the results obtained with our scheme LMe and the reference 
solution obtained by the explicit scheme at times t = 0.2, t = 0.4 and t = 0.8. Our 
scheme LMe is simulated with 50 and 200 grid points in space, whereas the explicit 
scheme is simulated with 2000 grid points. Again, we observe a good agreement 
between the two methods, in the kinetic regime. 

Example 4- Intermediate regime with non isotropic boundary conditions: 

a{v,w) = l, h{v) = v, fr{v) = Q, £ = 10-2. 

On Figures [7] and [8j we plot the results obtained with our scheme LMe, the schemes 
K, JPT, LMi with 10 gridpoints or 50 gridpoints in space, and compare them to 
the reference solution obtained with the explicit scheme using 10000 gridpoints in 
space, at time t = 0.4. It is clear that our scheme is in a good agreement with 



18 



M. LEMOU AND F. MEHATS 



the reference, even with coarse grids. There is no need in our scheme to discretize 
the boundary layer. Moreover, we recall that no artificial boundary condition is 
imposed, since our scheme is constructed in such a way that only the natural inflow 
boundary condition is used. We also see on Figures [7] and [8] that the numerical 
boundary value (at x = 0) accurately approximates the right value and we have 
the following remarkable property: the curve obtained by our scheme is close to 
the reference curve whatever the mesh size is. This property is fulfilled in all our 
numerical tests. 

Example 5. Diffusive regime with non isotropic boundary conditions: 

a{v,w) = l, fe{v) = v, Mv) = 0, £ = 10"^ 

On Figures [9] and [TOj we plot the results obtained with our scheme LMe, the schemes 
K, JPT, LMi with 50 gridpoints or 200 gridpoints in space, and compare them to 
the reference solution obtained using a scheme for the diffusion limiting equation 
(with the exact Chandrasekhar value at the left boundary), at time t = 0.4. Again, 
we observe that the results of our scheme LMe fit well with the reference curve. For 
the scheme K, we observe that the boundary value is very close to the right one, 
but the diffusive behavior is not correctly reproduced. For the schemes JPT and 
LMi, we observe that the uncorrectness of the diffusive regime is due to the fact 
that boundary condition is not the right one. 

On Figure [TT] we plot the density obtained by our scheme LMe using different 
numbers of space gridpoints: 10, 25, 50, 100 and 200, compared to the diffusion. 
We clearly observe that our scheme has a good behavior even with coarse grids. 

4.3. Collision operators with non constant cross sections. In this subsection, 
we provide some numerical simulations in cases where the cross section a is non 
constant. Various expressions of a{v,w) will be tested. In the first two examples. 



the cross section takes a special form allowing to derive a formula similar to (2.11) 
in order to compute numerically the associated generalized Chandrasekhar function. 
In these cases, we accurately know the exact value at the boundary in the diffusive 
regime, and this enables to compare our numerical results with this value. In the 
other two cases, we shall only compare our results with those obtained by the explicit 
scheme, in the intermediate regime e = 10-2. In all cases, we only consider a non 
isotropic boundary condition at x = 0: fe{v) = v, fr{v) = 0, which induces a 
boundary layer in the diffusive regime. 

Before presenting the numerical tests, let us give the generalized Chandrasekhar 
like formula corresponding to the specific cases where the cross section has the form 
a{v, w) = p{v)p{w), p being an even function. Following the method presented for 
instance in [H [13], one can derive the following formula. Let us denote the limiting 
boundary value at a; = for the diffusion problem by 

p,(0) = C K^{v)h{v)dv. 



Jo 

Then the kernel can be determined thanks to the generalized Chandrasekhar 
function as follows: 



K^iv) = -vH^{v) ^J^ piw)dwj i^J^ —dw ] , (4.1) 
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where H„ satisfies the nonhnear equation 



H^{v) = 1 + 



1 



w 



vp{w) + wp{v 



-H„{w)dw. (4.2) 



We shall compare numerically this kernel K(j{v) with its approximation Ka^{v) 
given by our scheme, which has been computed above and is given by (3.34). 



Example 6: a{v,w) = |f p/^|t(7p/^. 

We first plot on Figure [2] the approximated kernel Ka^^{v) given by (3.34), compared 
to the exact kernel Kf^{v) computed thanks to (|4.1[), (4.2). We observe that the 



kernel generated by our scheme is a very accurate approximation of the exact kernel. 




0.4 0.5 0.6 

velocity 



Figure 2. Example 6. Exact kernel Ka-{v) given by (4.2) and its 
approximation K^^^ defined by (3.34). 



We plot the density in the same three different regimes as above: e = 1 on Figure 



12 



e = 10~^ on Figure 13 and e = 10~^ on Figure 14 In all cases, the observations 
we made in the Subsection 14 . 2 1 for a constant collision kernel are still valid here and 
our curves are in good agreement with the reference. 



Example 7: a{v,w) = (1 — ^/^(l — 

We plot on Figure [3] the approximated kernel Ka^{v) given by (3.34), compared 
to the exact kernel K„[v) computed thanks to (4.1), (4.2). We observe again that 



the kernel generated by our scheme is a very accurate approximation of the exact 
kernel. 

We plot the density in the different regimes: e = 1 on Figure 15, e = 10^^ on 
Figure 16 and e = 10""^ on Figure 17 Again, our numerical tests provide a good 
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Figure 3. Example 7. Exact kernel K„{v) given by (4.2) and its 
approximation K'^^^ defined by (3.34) 



approximation of the exact result and we do not need any mesh refinement when a 
boundary layer appears. 

In the two following last examples, we consider cross sections for which formula 



analogous to (2.11) and (4.2) are not available and could be more complicated to 
obtain. We only test the intermediate regime e = 10~^, where the reference solution 
is computed thanks to the highly resolved explicit scheme. 



Example 8: a{v, w) = {v — wl^, intermediate regime e = 10~^. The density computed 
with our scheme is plotted on Figure 18 for 10 and 50 space gridpoints and is in 
good agreement with the reference computation. 



Example 9: a{v,w) = \v — w\~^'^, intermediate regime e = 10~^. The density 
computed with our scheme is plotted on Figure 19, again for 10 and 50 space 
gridpoints, and is also in good agreement with the reference computation. 



5. Conclusion 

In this paper, we have introduced a new strategy to develop the so-called Asymp- 
totic Preserving (AP) scheme for kinetic problems with possible boundary layers in 
the diffusion limit. This strategy is based on a new micro-macro decomposition. The 
macroscopic part of the distribution function is not the usual associated Maxwellian 
but a modified equilibrium which has the same incoming moments. This method is 
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then used to construct a numerical scheme which is AP inside the domain and pro- 
vides a very good approximation of the exact boundary condition, in both kinetic 
and diffusion regimes. This scheme only uses the natural inflow boundary condition, 
and in particular we do not solve the associated Milne problem and we do not need 
to inject the theoretical limiting boundary value given by the Chandrasekhar type 
formula (note moreover that such formula is not always available). 

We believe that our approach is sufficiently robust to be applied in other situa- 
tions. For example, we will explore this strategy in the case of hydrodynamic scaling 
with BGK type collision operators and then with more complex collision operators 
such as the Landau operator of plasma physics. Moreover, we emphasize that the 
approach can be extended naturally to multidimensional situations and this task 
will be achieved is a future work. 
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Figure 5. Example 2: comparison between the diffusion limit and 
our scheme LMe (50 and 200 space gridpoints), e = 10~^. Results 
at times t = 0.01, t = 0.1 and t = 0.4. 
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Figure 6. Example 3: comparison between the explicit scheme and 
our scheme LMe (50 and 200 space gridpoints), £ = 1. Results at 
times t = 0.2, t = 0.4 and t = 0.8. 




Figure 7. Example 4: comparison between the explicit scheme (10 
000 space gridpoints) and the schemes K, JPT, LMe and LMi (10 
space gridpoints each), e = 10~^. A zoom has been made on the 
axes: x G [0, 0.5]. 
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Figure 8. Example 4: comparison between the explicit scheme (10 
000 space gridpoints) and the schemes K, JPT, LMe and LMi (50 
space gridpoints each), e = 10~^. A zoom has been made on the 
axes: x G [0, 0.5]. 




Figure 9. Example 5: comparison between the diffusion (200 space 
gridpoints) and the schemes K, JPT, LMe and LMi (50 space grid- 
points each), £ = 10~^. A zoom has been made on the axes: x G 
[0,0.5]. 



M. LEMOU AND F. MEHATS 




0.45 



0.4 I ' ' ' ' ' ' ' ' ' 1 

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 

X 

Figure 10. Example 5: comparison between the diffusion (200 
space gridpoints) and tlie scliemes K, JPT, LMe and LMi (200 space 
gridpoints eacli), £ = 10"'*. A zoom has been made on the axes: 
X G [0,0.2]. 




Figure 11. Example 5: comparison between the diffusion (200 
space gridpoints) and our scheme LMe with 10, 25, 50, 100 and 
200 space gridpoints, £ = 10~^. A zoom has been made on the axes: 
X G [0,0.4]. 
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Figure 12. Example 6, kinetic regime e = 1: comparison between 
the explicit scheme (2000 space gridpoints) and our scheme LMe (50 
and 200 space gridpoints). Results at times t = 0.2, t = 0.4 and 
t = 0.8. 
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Figure 13. Example 6, intermediate regime e = 10~^: comparison 

between the explicit scheme (10 000 space gridpoints) and our scheme 
LMe (10 and 50 space gridpoints), at time t = 0.4. 
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Figure 14. Example 6, diffusive regime e = 10~^: comparison be- 
tween the diffusion and our scheme LMe (50 and 200 space grid- 
points), at time t = 0.4. 




Figure 15. Example 7, kinetic regime e = 1: comparison between 
the explicit scheme (2000 space gridpoints) and our scheme LMe (50 
and 200 space gridpoints). Results at times t = 0.2, t = 0.4 and 
t = 0.8. 
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Figure 17. Example 7, diffusive regime e = 10""^: comparison be- 
tween the diffusion and our scheme LMe (50 and 200 space grid- 
points), at time t = 0.4. 
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Figure 18. Example 8, intermediate regime e = 10~^: comparison 

between the explicit scheme (10 000 space gridpoints) and our scheme 
LMe (10 and 50 space gridpoints), at time t = 0.4. 
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Figure 19. Example 9, intermediate regime e = 10~^: comparison 
between the explicit scheme (10 000 space gridpoints) and our scheme 
LMe (10 and 50 space gridpoints), at time t = 0.4. 



